Eigen-vector approach for coil sensitivity maps estimation

ABSTRACT

A method for estimating a coil sensitivity map for a magnetic resonance (MR) image includes providing ( 61 ) a matrix A of sliding blocks of a 2D image of coil calibration data, calculating ( 62 ) a left singular matrix V ∥  from a singular value decomposition of A corresponding to τ leading singular values, calculating ( 63 ) P=V ∥ V ∥   H , calculating ( 64 ) a matrix S that is an inverse Fourier transform of a zero-padded matrix P, and solving ( 65 ) M H c r =(S r ) H c r  for c r , where c r  is a vector of coil sensitivity maps for all coils at spatial location r, and 
     
       
         
           
             
               M 
                
               
                 ( 
                 
                   
                     ( 
                     
                       
                         
                           1 
                         
                         
                           1 
                         
                         
                           … 
                         
                         
                           1 
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                       
                         
                           … 
                         
                         
                           … 
                         
                         
                           
                               
                           
                         
                         
                           … 
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                     
                     ) 
                   
                    
                   
                     ( 
                     
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                       
                         
                           1 
                         
                         
                           1 
                         
                         
                           … 
                         
                         
                           1 
                         
                       
                       
                         
                           … 
                         
                         
                           … 
                         
                         
                           
                               
                           
                         
                         
                           … 
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                     
                     ) 
                   
                    
                   
                       
                   
                    
                   … 
                    
                   
                       
                   
                    
                   
                     ( 
                     
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                       
                         
                           … 
                         
                         
                           … 
                         
                         
                           
                               
                           
                         
                         
                           … 
                         
                       
                       
                         
                           1 
                         
                         
                           1 
                         
                         
                           … 
                         
                         
                           1 
                         
                       
                     
                     ) 
                   
                 
                 ) 
               
             
             .

CROSS REFERENCE TO RELATED UNITED STATES APPLICATIONS

This application claims priority from “On the Eigen-Vector Approach forCoil Sensitivity Maps Estimation”, U.S. Provisional Application No.61/618,002 of Liu, et al., filed Mar. 30, 2012, “Coil Sensitivity mapsEstimation in Dynamic CMRI: Square Could be Better than Rectangular”,U.S. Provisional Application No. 61/724,023 of Liu, et al., filed Nov.8, 2012, and “An Eigen-Vector Approach for Coil Sensitivity Estimationin the 3D Scenario”, U.S. Provisional Application No. 61/724,001 of Liu,et al., filed Nov. 8, 2012, the contents of all of which are hereinincorporated by reference in their entireties.

TECHNICAL FIELD

This disclosure is directed to methods for estimating coil sensitivitymaps (CSM) of magnetic resonance imaging (MRI) apparatuses.

DISCUSSION OF THE RELATED ART

Parallel imaging makes use of multiple receiver coils to acquire theimage in parallel. It can be used to accelerate the image acquisition byexploiting the spatially varying sensitivities of the multiple receivercoils since each coil image is weighted differently by the coilsensitivity maps (CSM). The accurate estimation of coil sensitivity mapscan ensure the success of approaches that make use of the sensitivitymaps either in the reconstruction formulation or in coil combination.Some methods explicitly make use of the coil sensitivities inreconstruction, and the goodness of the CSM's is important to thequality of the reconstructed image. Other approaches implicitly make useof the coil sensitivities by performing an auto-calibrating coil-by-coilreconstruction, but the CSM's are needed if one wants to obtain the coilcombined complex image or the phase image.

The most common way to determine the sensitivity maps is to obtainlow-resolution pre-scans. However, when the object is not static, thesensitivity functions are different between pre-scan and under-sampledscans, and this could lead to reconstruction errors. To compensate forthis, joint estimation approaches have been proposed, however, theseapproaches usually have high computation cost and are restricted toexplicit reconstructions.

The eigenvector approach proposed tries to build a connection betweenimplicit and explicit approaches, by showing that the CSM can becomputed with the auto-calibrated coil-by-coil reconstruction. It hasbeen shown that the coil sensitivities can be computed as theeigenvector of a given matrix in the image space corresponding to theeigenvalue's ‘1”s. However, to the best of the inventor's knowledge: (1)the detailed mathematical derivations for the eigenvector approach arenot well understood; (2) the optimization criterion for computing theCSM is not very clear; and (3) there lacks an efficient approach.

SUMMARY

Exemplary embodiments of the invention as described herein generallyinclude methods for coil sensitivity maps (CSM) estimation for 2-D MRimages. Embodiments of the invention mathematically derive aneigenvector approach and show that the coil sensitivity maps at a givenspatial location can be obtained by solving a generalized eigenvaluesystem. Embodiments of the invention establish relationships between thegeneralized eigenvalue system and two related eigenvalue systems wherethe associated matrices are Hermitian. Embodiments of the invention canreduce the computational and storage costs and avoid the computation oflarge matrices by using equivalent representations. Experimental resultson simulated and real data sets show that, a method according to anembodiment of the invention can obtain CSM's with high quality and theresulting reconstructed images have fewer artifacts compared to thosereconstructed with CSM's obtained by a B1-Map.

According to an aspect of the invention, there is provided a method forestimating a coil sensitivity map for a magnetic resonance (MR) image,including providing a matrix A of sliding blocks of a 2D n_(x)×n_(y)image of coil calibration data, where

${A = \begin{pmatrix}a_{1,1} & a_{1,2} & \ldots & a_{1,n_{b}} \\a_{2,1} & a_{2,2} & \ldots & a_{2,n_{b}} \\\ldots & \ldots & \; & \ldots \\a_{n_{c,}1} & a_{n_{c},2} & \ldots & a_{n_{c},n_{b}}\end{pmatrix}},$

n_(c) is a number of coils, n_(b) is a number of sliding blocksextracted from the coil calibration data, and a_(i,j) is a k_(x)k_(y)×1column vector that represents a jth sliding block of an ith coil,determining a unit eigenvector α that maximizes

(S^(r))^(H)α,M^(H)α

, where S^(r) is an inverse Fourier transform of a zero-padded matrixP=V_(∥)V_(∥) ^(H) at spatial location r in the 2D n_(x)×n_(y) image,where

$V_{} = {\begin{bmatrix}{v_{1},} & {v_{2},\ldots \mspace{14mu},} & v_{\tau}\end{bmatrix} = \begin{pmatrix}v_{1,1} & v_{1,2} & \ldots & v_{1,\tau} \\v_{2,1} & v_{2,2} & \ldots & v_{2,\tau} \\\ldots & \ldots & \; & \ldots \\v_{n_{c},1} & v_{n_{c},2} & \ldots & v_{n_{c},\tau}\end{pmatrix}}$

is a matrix composed of left singular vectors of A corresponding to τleading singular values where v_(i,k) is a k_(x)k_(y)×1 column vectorthat is an i-th block in vector v_(k),

$M\left( {\begin{pmatrix}1 & 1 & \ldots & 1 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\begin{pmatrix}0 & 0 & \ldots & 0 \\1 & 1 & \ldots & 1 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\1 & 1 & \ldots & 1\end{pmatrix}} \right)$

is a n_(c)×k_(x)k_(y)n_(c) sparse matrix of the same size as S^(r), andfinding an optimizer c^(r) that maximizes a correlation between(S^(r))^(H)α and M^(H)α where the optimizer c^(r) is a column vector ofcoil sensitivities of a coil at spatial position r in the 2D n_(x)×n_(y)image.

According to another aspect of the invention, there is provided a methodfor estimating a coil sensitivity map for a magnetic resonance (MR)image, including providing a matrix A of sliding blocks of a 2Dn_(x)×n_(y) image of coil calibration data, where

${A = \begin{pmatrix}a_{1,1} & a_{1,2} & \ldots & a_{1,n_{b}} \\a_{2,1} & a_{2,2} & \ldots & a_{2,n_{b}} \\\ldots & \ldots & \; & \ldots \\a_{n_{c,}1} & a_{n_{c},2} & \ldots & a_{n_{c},n_{b}}\end{pmatrix}},$

n_(c) is a number of coils, n_(b) is a number of sliding blocksextracted from the coil calibration data, and a_(i,j) is a k_(x)k_(y)×1column vector that represents a jth sliding block of an ith coil,calculating a left singular matrix V_(∥) from a singular valuedecomposition of A, where A=VΣU^(H) and

$V_{} = {\begin{bmatrix}{v_{1},} & {v_{2},\ldots \mspace{14mu},} & v_{\tau}\end{bmatrix} = \begin{pmatrix}v_{1,1} & v_{1,2} & \ldots & v_{1,\tau} \\v_{2,1} & v_{2,2} & \ldots & v_{2,\tau} \\\ldots & \ldots & \; & \ldots \\v_{n_{c},1} & v_{n_{c},2} & \ldots & v_{n_{c},\tau}\end{pmatrix}}$

is a matrix of left singular vectors of A corresponding to τ leadingsingular values where v_(i,k) is a k_(x)k_(y)×1 column vector that is ani-th block in vector V_(k), calculating

$P = {{V_{}V_{}^{H}} = \left( {\begin{pmatrix}p_{1,1,1} & p_{1,1,2} & \ldots & p_{1,1,{k_{x}k_{y}}} \\p_{1,1,1} & p_{1,2,2} & \ldots & p_{1,2,{k_{x}k_{y}}} \\\ldots & \ldots & \ldots & \ldots \\p_{1,n_{c},1} & p_{1,n_{c},2} & \ldots & p_{1,n_{c},{k_{x}k_{y}}}\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}p_{n_{c},1,1} & \ldots & p_{n_{c},1,{k_{x}k_{y}}} \\p_{n_{c},2,1} & \ldots & p_{n_{c},2,{k_{x}k_{y}}} \\\ldots & \ldots & \ldots \\p_{n_{c},n_{c},1} & \ldots & p_{n_{c},n_{c},{k_{x}k_{y}}}\end{pmatrix}} \right)}$

where p_(i,j,t) is a k_(x)k_(y)×1 column vector, calculatingS_(i,j,t)=F^(H)(P_(t)(p_(i,j,t))) where F^(H) represents an inverseFourier transform and P_(t) represents a zero-padding operator, andsolving M^(H)c^(r)=(S^(r))^(H)c_(r) for c^(r), where c^(r) is a vectorof coil sensitivity maps for all coils at spatial location r,

$S^{r} = {\left( {\begin{pmatrix}s_{1,1,1}^{r} & s_{1,1,2}^{r} & \ldots & s_{1,1,{k_{x}k_{y}}}^{r} \\s_{1,2,1}^{r} & s_{1,2,2}^{r} & \; & s_{1,2,{k_{x}k_{y}}}^{r} \\\ldots & \ldots & \ldots & \ldots \\s_{1,n_{c},1}^{r} & s_{1,n_{c},2}^{r} & \ldots & s_{1,n_{c},{k_{x}k_{y}}}^{r}\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}s_{n_{c},1,1}^{r} & \ldots & s_{n_{c},1,{k_{x}k_{y\;}}}^{r} \\s_{n_{c},2,1}^{r} & \ldots & s_{n_{c},2,{k_{x}k_{y\;}}}^{r} \\\ldots & \ldots & \ldots \\s_{n_{c},n_{c},1}^{r} & \ldots & s_{n_{c},n_{c},{k_{x}k_{y\;}}}^{r}\end{pmatrix}} \right)\mspace{14mu} {and}}$$M = {\left( {\begin{pmatrix}1 & 1 & \ldots & 1 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\begin{pmatrix}0 & 0 & \ldots & 0 \\1 & 1 & \ldots & 1 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\1 & 1 & \ldots & 1\end{pmatrix}} \right).}$

According to a further aspect of the invention, the method includesfinding an optimizer c^(r) that maximizes a correlation between(S^(r))^(H)α and M^(H)α where the optimizer c^(r) is vector of coilsensitivity maps for all coils at spatial location r in the 2Dn_(x)×n_(y) image.

According to a further aspect of the invention, solvingM^(H)c^(r)=(S^(r))^(H)c^(r) comprises solving eigenvalue equations

${{\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}\beta} = {\overset{\_}{\lambda}\beta}},$

and (S_(r)(S^(r))^(H))γ={tilde over (λ)}γ, wherein λ and {tilde over(λ)} denote eigenvalues, β and γ denote eigenvectors.

According to a further aspect of the invention, the method includescalculating

${{S^{r}M^{H}\mspace{14mu} {using}\mspace{14mu} m_{i,j}} = {{\sum\limits_{t = 1}^{k_{x}k_{y}}\; s_{i,j,t}} = {{\sum\limits_{t = 1}^{k_{x}k_{y}}\; {F^{H}\left( {P_{t}\left( p_{i,j,t} \right)} \right)}} = {F^{H}\left( {\sum\limits_{t = 1}^{k_{x}k_{y}}\; {P_{t}\left( p_{i,j,t} \right)}} \right)}}}},$

where m_(i,j) is an n_(x)n_(y)×1 column vector containing the entries ofm_(i,j) ^(r) at all spatial locations, and m_(i,j) ^(r) is the i,jelements of the 2D MR image at spatial location r.

According to a further aspect of the invention, the method includescalculating S^(r)(S^(r))^(H) asS^(r)(S^(r))^(H)=k_(x)k_(y)Z^(r)(Z^(r))^(H), where

$Z^{r} = \begin{pmatrix}z_{1,1}^{r} & z_{1,2}^{r} & \; & z_{1,\tau}^{r} \\z_{2,1}^{r} & z_{2,2}^{r} & \; & z_{2,\tau}^{r} \\\; & \; & \; & \; \\z_{n_{c},1}^{r} & z_{n_{c},2}^{r} & \; & z_{n_{c},\tau}^{r}\end{pmatrix}$

is an n_(c)×τ matrix with the j-row and k-th column being Z_(j,k) ^(r),Z_(j,k) ^(r) denotes the r-th entry of Z_(j,k)=F^(H) (P_(c)(v_(j,k)))and P_(c) denotes the zero-padding operator that maps a center ofv_(j,k) to a center of z_(j,k).

According to a further aspect of the invention, the method includescomputing the i-th row and j-th column of Z^(r)(Z^(r))^(H) as

$\sum\limits_{k = 1}^{\tau}\; {z_{i,k}^{r}\overset{\_}{z_{i,k}^{r}}}$

where z_(i,k) ^(r) denotes the i-th row and k-th column in Z^(r).

According to a further aspect of the invention, the coil calibrationdata is acquired from a center square of frequency space data.

According to a further aspect of the invention, the matrix A isconstructed from calibration data of size c_(x)×c_(y)×c_(z) by gatheringsliding blocks of kernel size k_(x)×k_(y)×k_(z), where A has[(c_(x)−k_(x)+1)×(c_(y)−k_(y)+1)×(c_(z)−k_(z)+1)] columns andk_(x)k_(y)k_(z)n_(c) rows.

According to another aspect of the invention, there is provided anon-transitory program storage device readable by a computer, tangiblyembodying a program of instructions executed by the computer to performthe method steps for estimating a coil sensitivity map for a magneticresonance (MR) image.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a table summarizing the notation used in this disclosure,according to an embodiment of the invention.

FIG. 2 illustrates zero padding for the case k_(x)=k_(y)=3, according toan embodiment of the invention.

FIG. 3 illustrates an eigenvector approach where each k-space point isrepresented as a linear combination of the k-space points of all thecoils in its neighborhood, according to an embodiment of the invention.

FIG. 4 illustrates the principle of an eigenvector approach of EQ. (13),according to an embodiment of the invention.

FIG. 5 is a table summarizing the value of mean absolute correlation(mac) under different number of kept singular vectors and differentkernel sizes, according to an embodiment of the invention.

FIG. 6 is a flowchart of an eigenvector approach to estimating coilsensitivity maps, according to an embodiment of the invention.

FIGS. 7( a)-(d) depict CSM's obtained by various approaches, accordingto an embodiment of the invention.

FIGS. 8( a)-(c) are graphs of the mean absolute correlation as afunction of different number of kept vectors and different kernel sizes,according to an embodiment of the invention.

FIG. 9 illustrates a reconstruction of a 2D cardiac image, according toan embodiment of the invention.

FIGS. 10( a)-(b) illustrates a reconstruction of a 3D abdomen image,according to an embodiment of the invention.

FIG. 11 illustrates a typical Cartesian sampling pattern used in dynamiccardiac MR imaging, according to an embodiment of the invention.

FIGS. 12( a)-(b) illustrate the influence of the amount of usedcalibration data on the reconstruction, according to an embodiment ofthe invention.

FIG. 13 compares the square and the rectangular calibration data interms of the mac criterion, according to an embodiment of the invention.

FIG. 14 illustrates the conversion of the calibration data into acalibration matrix for 3D CSM estimation, according to an embodiment ofthe invention.

FIGS. 15( a)-(h) illustrate a comparison of a 2D eigenvector approachand a 3D eigenvector approach according to an embodiment of theinvention.

FIG. 16 is a block diagram of an exemplary computer system forimplementing an eigenvector approach to coil sensitivity maps estimationfor 2-D MR images, according to an embodiment of the invention.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

Exemplary embodiments of the invention as described herein generallyinclude systems and methods for an eigenvector approach to coilsensitivity maps (CSM) estimation for 2-D MR images. Accordingly, whilethe invention is susceptible to various modifications and alternativeforms, specific embodiments thereof are shown by way of example in thedrawings and will herein be described in detail. It should beunderstood, however, that there is no intent to limit the invention tothe particular forms disclosed, but on the contrary, the invention is tocover all modifications, equivalents, and alternatives falling withinthe spirit and scope of the invention.

As used herein, the term “image” refers to multi-dimensional datacomposed of discrete image elements (e.g., pixels for 2-dimensionalimages and voxels for 3-dimensional images). The image may be, forexample, a medical image of a subject collected by computer tomography,magnetic resonance imaging, ultrasound, or any other medical imagingsystem known to one of skill in the art. The image may also be providedfrom non-medical contexts, such as, for example, remote sensing systems,electron microscopy, etc. Although an image can be thought of as afunction from R³ to R or R⁷, the methods of the inventions are notlimited to such images, and can be applied to images of any dimension,e.g., a 2-dimensional picture or a 3-dimensional volume. For a 2- or3-dimensional image, the domain of the image is typically a 2- or3-dimensional rectangular array, wherein each pixel or voxel can beaddressed with reference to a set of 2 or 3 mutually orthogonal axes.The terms “digital” and “digitized” as used herein will refer to imagesor volumes, as appropriate, in a digital or digitized format acquiredvia a digital acquisition system or via conversion from an analog image.

Notation:

Throughout this disclosure, scalars are denoted by italic lowercaseletters (e.g., n_(c), n_(x), n_(y), k_(x), k_(y)), vectors by boldfacelowercase letters (e.g., x_(i), c_(i)), and matrices by italic uppercaseletters (e.g., S S^(r)). The main notations used in this paper aresummarized in Table I, shown in FIG. 1.

Let { }^(H) denote the complex-conjugate transpose operator, x denotethe complex-conjugate of a scalar x, and x be the complex-conjugate of avector x.

Let i and j denote the indices for the coils in the interval [1, n_(c)],r denote the spatial location in the 2D image of size n_(x)×n_(y), and tdenote the spatial location in the image block of size k_(x)×k_(y).

For an n_(x)n_(y)×1 column vector x_(i), representing an image byconcatenating its columns, x_(i) ^(r) corresponds to its value at thespatial location r, where r goes from 1 to n_(x)n_(y). For ak_(x)k_(y)×1 column vector v_(i,k), representing an image block byconcatenating its columns, v_(i,k) ^(t) denotes its value at the spatiallocation t, where t goes from 1 to k_(x)k_(y).

Let y=P_(t)(x) denote the zero padding operator defined as follows. Theinput x, a column vector of size k_(x)k_(y)×1, represents a k_(x)×k_(y)image block by concatenating its columns; t, a scalar in the interval[1, k_(x)k_(y)], determines which value in x to be put at the center ofthe matrix of size n_(x)×n_(y); the output y is an n_(x)×n_(y)×1 columnvector representing an n_(x)×n_(y) matrix, and its non-zero entries arethe values of x. In particular, we let y=P_(c)(x) denote the zeropadding operator that put the center point of x to the center locationof y. For an image of size n_(x)×n_(y), define the index of its centerlocation as

${\left\lceil \frac{n_{x} + 1}{2} \right\rceil + {\left( {\left\lceil \frac{n_{y} + 1}{2} \right\rceil - 1} \right)n_{x}}};$

and for an image block of size k_(x)×k_(y), define the index of itscenter location as

$\left\lceil \frac{k_{x} + 1}{2} \right\rceil + {\left( {\left\lceil \frac{k_{y} + 1}{2} \right\rceil - 1} \right){k_{x}.}}$

Embodiments of the present disclosure assume that n_(x) and n_(y) areeven; but k_(x) and ky can be either even or odd.

FIG. 2 illustrates zero padding for the case k_(x)=k_(y)=3. The ninesquares correspond to the input x, the square 21 is placed at the centerof the n_(x)×n_(y) image, t in

determines which value in x to be put at the center of the 2D image; theremaining entries are filled with zero except for the nine squares; andthe three rectangles 22, 23, 24 denote the outputs of P₁, P₅, and P₉,respectively.

Eigenvector Approach

One eigenvector approach is illustrated in FIG. 3, where each k-spacepoint is represented as a linear combination of the k-space points ofall the coils in its neighborhood. Let x_(i) ^(r) denote the value ofx_(i) at the spatial location r. Let w_(i,j) be the convolution kernelfor representing the k-space point of the i-th coil by the j-th coil.Mathematically, x_(i) ^(r) is estimated as:

$\begin{matrix}{{x_{i}^{r} = {\sum\limits_{j = 1}^{n_{c}}\; {w_{i,j}*x_{j}}}},{\forall r},} & (1)\end{matrix}$

where * denotes the convolution operator. Define

q _(i,j) =F ^(H)(P _(C)(w _(i,j))).  (2)

q_(i,j) is the result of applying an inverse Fourier transformation tothe convolution kernel w_(i,j) with zero padding, illustrated in FIG. 2.Applying the inverse Fourier transform to both sides of the equation inFIG. 3, one has

$\begin{matrix}{{F^{H}\left( x_{i} \right)} = {\sum\limits_{j = 1}^{n_{c}}\; {q_{ij} \otimes {F^{H}\left( x_{j} \right)}}}} & (3)\end{matrix}$

where the convolution theorem has been applied. F^(H)(x_(i))=m

c_(i) can be inserted into EQ. (3) to obtain

$\begin{matrix}{{c_{i} = {\sum\limits_{j = 1}^{n_{c}}\; {q_{ij} \otimes c_{j}}}},\mspace{14mu} {i = 1},2,\ldots \mspace{14mu},{n_{c}.}} & (4)\end{matrix}$

EQ. (4) indicates that all the entries in this equation can be decoupledwith respect to the spatial location. Let c_(i) ^(r), q_(i,j) ^(r) andc_(j) ^(r) denote the entries of c_(i), q_(i,j), and c_(j) at thespatial location r, respectively. EQ. (4) can be rewritten as

$\begin{matrix}{{c_{i}^{r} = {\sum\limits_{j = 1}^{n_{c}}\; {q_{i,j}^{r}c_{j}^{r}}}},\mspace{14mu} {i = 1},2,\ldots \mspace{14mu},n_{c},{\forall{r.}}} & (5)\end{matrix}$

Denote Q^(r) as a matrix with the value at the i-row and the j-th columnbeing q_(i,j) ^(r). Let c^(r) be an n_(c)×1 column vector composed ofc_(i)'s at the spatial location r. EQ. (5) can be rewritten as thefollowing compact matrix format:

c ^(r) =Q ^(r) c ^(r) ,∀r  (6)

FIG. 3 illustrates the principle of CSM: the three rectangles denote thek-space data of the i-th coil, the first coil, and the last coil,respectively; x_(i) ^(r) the square 31 in the left plot, is estimated asthe linear combination of the k-space points of all the coils in itsneighborhood of size k_(x)×k_(y) (here k_(x)=k_(y)=3); the nine smallsquares 32 in X₁ and 33 in X_(nc) correspond to the values in theconvolution kernel w_(i,l) and w_(i,n) _(c) .

In an eigenvector approach according to an embodiment of the invention,sliding blocks of size k_(x)×k_(y)) are gathered from the calibrationdata to form a matrix:

$\begin{matrix}{A = \begin{pmatrix}a_{1,1} & a_{1,2} & \ldots & a_{1,n_{b}} \\a_{2,1} & a_{2,2} & \ldots & a_{2,n_{b}} \\\ldots & \ldots & \; & \ldots \\a_{n_{c},1} & a_{n_{c},2} & \ldots & a_{n_{c},n_{b}}\end{pmatrix}} & (7)\end{matrix}$

where n_(b) denotes the number of sliding blocks extracted from thecalibration data, and a_(i,k), which is a k_(x)k_(y)×1 column vector,denotes the k-th sliding block of the i-th coil.

Let

A=VΣU ^(H)  (8)

be the Singular Value Decomposition of A. Denote

$\begin{matrix}{V_{} = {\left\lbrack {v_{1},v_{2},\ldots \mspace{14mu},v_{\tau}} \right\rbrack = \begin{pmatrix}v_{1,1} & v_{1,2} & \ldots & v_{1,\tau} \\v_{2,1} & v_{2,2} & \ldots & v_{2,\tau} \\\ldots & \ldots & \; & \ldots \\v_{n_{c},1} & v_{n_{c},2} & \ldots & v_{n_{c},\tau}\end{pmatrix}}} & (9)\end{matrix}$

as a matrix composed of the left singular vectors of A corresponding tothe leading τ singular values. Here, v_(i,k), a k_(x)k_(y)×1 columnvector, is the i-th block in the vector v_(k). It can be observed fromEQS. (7) and (9) that, A and V_(∥) share a similar structure in that:(1) the block a_(i,k) and v_(i,k) have the same size (k_(x)k_(y)×1); and(2) the number of such blocks for each column is same.

Let P=V_(∥)V_(∥) ^(H), and denote

$\begin{matrix}{P = \left( {\begin{pmatrix}p_{1,1,1} & p_{1,1,2} & \ldots & p_{1,1,{k_{x}k_{y}}} \\p_{1,1,1} & p_{1,2,2} & {\; \ldots} & p_{1,2,{k_{x}k_{y}}} \\\ldots & \ldots & \ldots & \ldots \\p_{1,n_{c},1} & p_{1,n_{c},2} & \ldots & p_{1,n_{c},{k_{x}k_{y}}}\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}p_{n_{c},1,1} & \ldots & p_{n_{c},1,{k_{x}k_{y\;}}} \\p_{n_{c},2,1} & \ldots & p_{n_{c},2,{k_{x}k_{y\;}}} \\\ldots & \ldots & \ldots \\p_{n_{c},n_{c},1} & \ldots & p_{n_{c},n_{c},{k_{x}k_{y\;}}}\end{pmatrix}} \right)} & (10)\end{matrix}$

where p_(i,j,t) is a k_(x)k_(y)×1 column vector. Embodiments have used asimilar block representation for P as A and V_(∥) in that they all haven_(c) “rows”. Note that the right hand side of EQ. (10) also representsP^(H), as P is Hermitian (i.e., P=P^(H)).

As V_(∥) and V have orthogonal column vectors,

V _(∥) V _(∥) ^(H) A=V _(∥) V _(∥) ^(H) ΣU ^(H) =[V∥0]ΣU ^(H) ≈A,  (11)

where I denotes an identity matrix. When τ equals the rank of A, the “≈”becomes “=”. EQ. (11) indicates that, each column of A lies in the spacespanned by orthogonal projection V_(∥)V_(∥) ^(H) when τ is set to therank of A.

Next, embodiments generate EQ. (11) for all k-space blocks.Specifically, let x_(i,r) be a k_(x)k_(y)×1 column vector that denotesthe k-space block of the i-th coil at the spatial location r. Assumingthat the k-space blocks lie in the space spanned by V_(∥),

$\begin{matrix}{{\begin{pmatrix}x_{1,r} \\x_{2,r} \\\ldots \\x_{n_{c},r}\end{pmatrix} = {V_{}{V_{}^{H}\begin{pmatrix}x_{1,r} \\x_{2,r} \\\ldots \\x_{n_{c},r}\end{pmatrix}}}},{\forall{r.}}} & (12)\end{matrix}$

The case when EQ. (12) holds approximately will be discussed in the nextsubsection.

Picking a column of P defined in EQ. (10) and making use ofP=P^(H)=V_(∥)V_(∥) ^(H) P, EQ. (12) yields

$\begin{matrix}{{x_{i,r}^{t} = {\sum\limits_{j = 1}^{n_{c}}\; {\left( p_{i,j,t} \right)^{H}x_{i,r}}}},{\forall r},} & (13)\end{matrix}$

where X_(i,r) ^(t) denotes the t-th element in the vector x_(i,r). FIG.4 illustrates the principle of an eigenvector approach according to anembodiment of the invention of EQ. (13) with k_(x)=k_(y)=3: x_(i,r) ^(t)denoted as the square 41 in the left plot, is represented by thesummation of the convolutions between the k-space data of the j-th coilwith kernels (p_(i:j:t))^(H) for j=1, 2, . . . , n_(c); the plots in thefirst, second, and third row correspond to t=1, 5, 9, respectively; andthe nine small squares 42 in X₁ and 43 in X_(nc), which denote(p_(i,j,t))^(H), constitute the convolution kernel. Note that, by movingthe sliding block x_(i,r) over all the spatial locations r=1, 2, . . . ,n_(x)n_(y), the left hand side of EQ. (13) covers the k-space data ofthe i-th coil; and the right hand side of EQ. (13) corresponds to thesummation of the convolutions between the k-space data of the j-th coiland the kernel (p_(i:j:t))^(H) for j=1, . . . , n_(c).

For a given i and t, applying the inverse Fourier transformation to bothsides of EQ. (13) yields

$\begin{matrix}{{c_{i} = {\sum\limits_{j = 1}^{n_{c}}\; {\overset{\_}{s_{i,j,t}} \otimes c_{j}}}},} & (14)\end{matrix}$

where i=1, 2, . . . , n_(c); t=1, 2, . . . , k_(x)k_(y), and

s _(i,j,t) =F ^(H)(P _(t)(p _(i,j,t)))  (15)

is the result of applying the inverse Fourier transform to p_(i,j,t)with proper zero padding, and s_(i,j,t) denotes the complex conjugate ofs_(i,j,t).

EQ. (14) indicates that the values at different spatial locations can bedecoupled. Let c_(i) ^(r), s_(i,j,t) ^(r) and c_(j) ^(r) denote theentries of c_(i), s_(i,j,t), and c_(j) at the spatial location r,respectively. Embodiments can rewrite EQ. (14) as

$\begin{matrix}{{c_{i}^{r} = {\sum\limits_{j = 1}^{n_{c}}\; {\overset{\_}{s_{i,j,t}^{r}}c_{j}^{r}}}},} & (16)\end{matrix}$

where i=1, 2, . . . , n_(c); t=1, 2, . . . , k_(x)k_(y).

Let c^(r) be a column vector composed of c_(i)'s at the spatial locationr. Then,

$\begin{matrix}{{S^{r} = \left( {\begin{pmatrix}s_{1,1,1}^{r} & s_{1,1,2}^{r} & \cdots & s_{1,1,{k_{x}k_{y}}}^{r} \\s_{1,2,1}^{r} & s_{1,2,2}^{r} & \; & s_{1,2,{k_{x}k_{y}}}^{r} \\\cdots & \cdots & \cdots & \cdots \\s_{1,n_{c},1}^{r} & s_{1,n_{c},2}^{r} & \cdots & s_{1,n_{c},{k_{x}k_{y}}}^{r}\end{pmatrix}{\cdots \begin{pmatrix}s_{n_{c},1,1}^{r} & \cdots & s_{n_{c},1,{k_{x}k_{y}}}^{r} \\s_{n_{c},2,1}^{r} & \cdots & s_{n_{c},2,{k_{x}k_{y}}}^{2} \\\cdots & \cdots & \cdots \\s_{n_{c},n_{c},1}^{r} & \cdots & s_{n_{c},n_{c},{k_{x}k_{y}}}^{r}\end{pmatrix}}} \right)}{Let}} & (17) \\{M = \left( {\begin{pmatrix}1 & 1 & \cdots & 1 \\0 & 0 & \cdots & 0 \\\cdots & \cdots & \; & \cdots \\0 & 0 & \cdots & 0\end{pmatrix}\begin{pmatrix}0 & 0 & \cdots & 0 \\1 & 1 & \cdots & 1 \\\cdots & \cdots & \; & \cdots \\0 & 0 & \cdots & 0\end{pmatrix}{\cdots \begin{pmatrix}0 & 0 & \cdots & 0 \\0 & 0 & \cdots & 0 \\\cdots & \cdots & \; & \cdots \\1 & 1 & \cdots & 1\end{pmatrix}}} \right)} & (18)\end{matrix}$

be a sparse matrix of the same size as S^(r) (i.e.,n_(c)×k_(x)k_(y)n_(c)). It can be verified that MM^(H)=k_(x)k_(y)I,where I denotes an identity matrix. EQ. (16) can be rewritten as

M ^(H) c ^(r)=(S ^(r))^(H) c ^(r).  (19).

EQ. (19) indicates that c^(r) is the eigenvector of a generalizedeigenvalue equation

(S ^(r))^(H) α=λM ^(H)α,  (20)

corresponding to the eigenvalue 1, where α is an n_(c)×1 column vectorrepresenting the generalized eigenvector and λ is a scalar representingthe generalized eigenvalue.

Since S^(r) and M are n_(c)×k_(x)k_(y)n_(c) non-square matrices, it isunclear (1) whether there exist a real eigenvalue λ and (complex)eigenvector α that satisfy EQS. (19) and (20); and (2) how one can solveEQ. (19) efficiently and robustly. The above situation is furthercomplicated by the fact that EQ. (12) usually only holds approximatelyas τ is usually set less than or equal to the rank of A, even for thecalibration data. Therefore, the left and right sides of EQS. (13),(14), (16), and (19) may only be approximately equal

To deal with the aforementioned issues, embodiments of the inventionsolve the following two eigenvalue equations:

$\begin{matrix}{{{\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}\beta} = {\overset{\_}{\lambda}\beta}},} & (21) \\{{{\left( {S^{r}\left( S^{r} \right)}^{H} \right)\gamma} = {\overset{\sim}{\lambda}\gamma}},} & (22)\end{matrix}$

where λ and λ denote the eigenvalues, β and γ denote the eigenvectors.Specifically, embodiments of the invention can compute c^(r) as theeigenvector of EQ. (21) corresponding to the largest eigenvalue. Inother words, c^(r) is the solution to the following equation

$\begin{matrix}{c^{r} = {\underset{{{\alpha \text{:}}\mathop{\text{||}}\alpha ||_{2}} = 1}{\arg \mspace{14mu} \max}{\langle{{\left( S^{r} \right)^{H}\alpha},{M^{H}\alpha}}\rangle}}} & (23)\end{matrix}$

and c^(r) is the optimizer that maximizes the correlation between(S^(r))^(H)α and M^(H)α, i.e., the correlation between the left andright hand sides of EQ. (19). Embodiments of the invention may refer tothis approach as MACO, which stands for maximal correlation. Inaddition, EQ. (23) also applies to the case in which EQ. (19) holdsexactly.

It will now be shown that M(S^(r))^(H) is Hermitian:

Theorem 1: The matrix M(S^(r))^(H) is Hermitian, i.e.,M(S^(r))^(H)=S^(r)M^(H). Thus, the eigenvalues of EQ. (21) are all real.

The proof of Theorem 1 is provided in the Appendix, below. Theeigen-decomposition of

$\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}$

has the following features: (1) the matrix is Hermitian, and theexistence of real eigenvalues is guaranteed; (2) there existsnumerically stable and efficient approaches for the eigenvaluedecomposition of a Hermitian matrix; and (3) the size of M(S^(r))^(H),n_(c)×n_(c), is much smaller than S^(r). In addition, experimentssuggest that M(S^(r))^(H) may be positive semi-definite.

Next, the relationships among EQS (20), (21), and (22) will be defined.

Theorem 2: For any λ and α satisfying the generalized eigenvalueequation of EQ. (20), λ is real, and

$\begin{matrix}{{{\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}\alpha} = {\lambda\alpha}},} & (24)\end{matrix}$

i.e., λ and α are also the eigenvalue and eigenvector of EQ, (21),respectively.

The proof of Theorem 2 is given in the Appendix, below. Theorem 2 showsthat any λ and α that satisfy EQ. (19) can be computed by applying theeigenvalue decomposition to the square Hermitian matrix

$\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}.$

Theorem 3: For any λ and α satisfying the generalized eigenvalueequation of EQ. (20), then

$\begin{matrix}{{{\frac{{S^{r}\left( S^{r} \right)}^{H}}{k_{x}k_{y}}\alpha} = {\lambda^{2}\alpha}},} & (25)\end{matrix}$

i.e., λ² and α are also the eigenvalue and eigenvector of EQ. (22),respectively. In addition, λ² and α are the eigenvalue and eigenvectorof EQ. (25), if and only if λ² and α are the eigenvalue and eigenvectorof

$\begin{matrix}{{{\left( \frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}} \right)^{2}\alpha} = {\lambda^{2}\alpha}},} & (26)\end{matrix}$

The proof of Theorem 3 is given in the Appendix, below. This theoremshows that: (1) if a generalized eigenvalue equation has solutions, itcan also be solved by EQ. (22); (2) the eigenvalue decomposition of

$\frac{{S^{r}\left( S^{r} \right)}^{H}}{k_{x}k_{y}}$

has the same eigenvalues and eigenvectors as

$\left( \frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}} \right)^{2};$

and (3) the eigenvector decomposition of

$\frac{{S^{r}\left( S^{r} \right)}^{H}}{k_{x}k_{y}}$

is also the eigenvector of

$\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}},$

as

$\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}$

is Hermitian. If the conjecture that

$\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}$

is positive semidefinite is true, then

$\frac{{S^{r}\left( S^{r} \right)}^{H}}{k_{x}k_{y}}$

has the same eigenvalues as

$\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}.$

A naive implementation would be to compute S^(r) for all the spatiallocations and then apply eigenvalue decomposition to EQ. (21) to obtainc^(r). However, the computational cost is n_(c) ²k_(x)k_(y) 2D Fouriertransformations, leading to a computation cost of O(n_(e)²k_(x)k_(y)n_(x)n_(y) log(n_(x)n_(y))) and a storage cost of O(n_(c)²k_(x)k_(y)n_(x)n_(y)). This can be too demanding for both computationand storage. For example, when n_(x)=n_(y)=256, k_(x)=k_(y)=6, andn_(c)=32, it costs about 10 GB for storage using single precision inMatlab and O(10⁹) floating operations for computation. When the sizekernel size k_(y)×k₃, is large, the computational and storage costs areeven higher.

To deal with the above issues, a MACO approach according to embodimentsof the invention need only to compute S^(r)M^(H) or S^(r)(S^(r))^(H). Auseful operation for efficiently computing S^(r)M^(H) is

$\begin{matrix}{{m_{i,j}^{r} = {\sum\limits_{t = 1}^{k_{x}k_{y}}s_{i,j,t}^{r}}},{\forall r},} & (27)\end{matrix}$

where m_(i,j) ^(r) is the i,j elements of the 2D MR image at spatiallocation r, which can be rewritten as

$\begin{matrix}{{m_{i,j} = {\sum\limits_{t = 1}^{k_{x}k_{y}}s_{i,j,t}}},} & (28)\end{matrix}$

where m_(i,j) is an n_(x)n_(y)×1 column vector containing the entries ofm_(i,j) ^(r) at all spatial locations. Incorporating EQ. (15):

$\begin{matrix}{{m_{i,j} = {{\sum\limits_{t = 1}^{k_{x}k_{y}}{F^{H}\left( {P_{t}\left( p_{i,j,t} \right)} \right)}} = {F^{H}\left( {\sum\limits_{t = 1}^{k_{x}k_{y}}{P_{t}\left( p_{i,j,t} \right)}} \right)}}},} & (29)\end{matrix}$

which shows that only one Fourier transformation is needed to obtainn_(i,j). As i and j goes from 1 to n_(c) and S^(r)M^(H) is Hermitian,

$\frac{n_{c}\left( {n_{c} + 1} \right)}{2}$

Fourier transformation are needed. With m_(i,j), it can be shown thatthe i-th row and j-th column of the matrix S^(r)M^(H) is m_(i,j) ^(r).The computation cost is O(n_(c) ²n_(x)n_(y) log(n_(x)n_(y))), and thememory cost is O(n_(e) ²n_(x)n_(y)).

When τ is small, EQ (22) can be more efficiently computed than EQ. (21).Let

z _(j,k) =F ^(H)(P _(c)(v _(j,k)))  (30)

and Z_(j,k) ^(r) denotes the r-th entry of z_(i,k). Let Z^(r) be ann_(c)×τ matrix with the j-row and k-th column being z_(j,k) ^(r), i.e.,

$\begin{matrix}{Z^{r} = \begin{pmatrix}z_{1,1}^{r} & z_{1,1}^{r} & \; & z_{1,\tau}^{r} \\z_{2,1}^{r} & z_{2,2}^{r} & \; & z_{2,\tau}^{r} \\\; & \; & \; & \; \\z_{n_{c},1}^{r} & z_{n_{c},2}^{r} & \; & z_{n_{c},\tau}^{r}\end{pmatrix}} & (31)\end{matrix}$

Theorem 4: S^(r)(S^(r))^(H) can be computed as:

S ^(r)(S ^(r))^(H) =k _(x) k _(y) Z ^(r)(Z ^(r))^(H)  (32)

The proof of Theorem 4 can be found in the Appendix, below.

By transforming the computation of S^(r)(S^(r))^(H) to Z^(r)(Z^(r))^(H)in EQ. (32), the number of 2-D Fourier transforms can be reduced fromn_(c) ²k_(x)k_(y) n_(c)τ. As τ is typically far less thann_(c)k_(x)k_(y), this can save the computation cost. In addition, thestorage cost can be reduced from O(n_(c) ²k_(x)k_(y)n_(x)n_(y)) toO(n_(c) ²k_(x)k_(y)), when Z^(r)(Z^(r))^(H) is implemented as follows.An implementation according to an embodiment of the invention is basedon the observation that the i-th row and j-th column of Z^(r)(Z_(r))^(H)can be computed as

$\begin{matrix}{\sum\limits_{k = 1}^{\tau}{z_{i,k}^{r}\overset{\_}{z_{i,k}^{r}}}} & (33)\end{matrix}$

where z_(i,k) ^(r) denotes the i-th row and k-th column in Z^(r). Thisshows that, instead of pre-computing Z^(r), embodiments of the inventioncan compute one column of Z^(r) each time.

Table II, shown in FIG. 5, summarizes the computation cost for computingS^(r), M(S^(r))^(H), and S^(r)(S^(r))^(H). It can be seen that thecomputational and storage costs of computing M(S^(r))^(H) andS^(r)(S^(r))^(H) are much slower than for S^(r). In addition, when

${\tau < \frac{n_{c} + 1}{2}},$

the computation cost of S^(r)(S^(r))^(H) is lower than M(S^(r))^(H).

FIG. 6 is a flowchart of an eigenvector approach to estimating coilsensitivity maps, according to an embodiment of the invention. Referringnow to the figure, an eigenvector approach begins at step 61 byproviding the matrix A of sliding blocks of a 2D n_(x)×n_(y) image ofcoil calibration data as defined in EQ. 7. At step 62, the left singularmatrix V_(∥) of EQ. (9) is calculated from the singular valuedecomposition of A of EQ. (8), and matrix P is calculated at step 63from V_(∥) using P=V_(∥)V_(∥) ^(H). At step 64, the elements s_(i,j,t)^(r) of matrix S^(r) are calculated by applying an inverse Fouriertransform to a zero-padded P, and at step 65, the eigenvalue equationM^(H)c^(r)=(S^(r))^(H)c^(r) of EQ. (19) is solved using the eigenvaluesystems of EQS. (21) and (22). The system of EQ. (21) can be solvedusing EQS. (27), (28), and (29) at step 66, and the system of EQ. (22)can be solved using EQ. (32), at step 67.

EXPERIMENTS

Simulation:

Experiments were performed on a simulated phantom. A 2D brain phantomand coil sensitivity maps of eight coil channels, were generated usingpublicly available codes, and the size of the brain phantom is set ton_(x)=n_(y)=256. The simulated eight coil sensitivity maps are shown inFIG. 7( a). With the simulated brain phantom and the simulatedsensitivity maps, k-space data was generated for the eight coils. Toestimate the coil sensitivity maps, the center 32×32 k-space data wasused as the calibration data. The estimated CSM's of a MACO approachaccording to an embodiment of the invention are shown in FIG. 7( b),where the kernel size is set to k_(x)=k_(y)=6 and τ=55. For comparison,FIG. 7( c) depicts CSM's estimated by the B1-Map approach with spatialscaling, and FIG. 7( d) depicts CSM's obtained by dividing the lowresolution coil images by the sum-square images. It can be observedthat, a MACO approach according to an embodiment of the inventiongenerates CSM's that are very close to the ground truth. Further, theestimated sensitivity maps in FIG. 7( d) in fact correspond to: (1)those estimated by the B1-Map approach if the spatial smoothing is notapplied; and (2) those estimated by a MACO approach according to anembodiment of the invention when the calibration kernel size is set tothe size of the calibration data.

Next, CSM's estimated by a MACO approach according to an embodiment ofthe invention were quantitatively compared with the ground truth and theB1-Map approach. Let c*^(r) be the ground truth coil sensitivities forall the coils at the spatial location r, assuming that c*^(r) has aunique length, i.e., ∥c*^(r)∥=1. Note that, in both a MACO approachaccording to an embodiment of the invention and a B1-Map, if c^(r) isthe obtained eigenvector, then ηc^(r) is also a valid eigenvector forthe corresponding eigenvalue if η is a complex scalar with unit length(i.e., |η|=1). The similarity between the estimated c^(r) and the groundtruth c*^(r) was measured using the absolute correlation defined as:

ρ(c ^(r) ,*c ^(r))=

c* ^(r) ,c ^(r)

|,  (34)

As ∥c*^(r)∥=∥c^(r)∥=1, if ρ(c^(r),c*^(r))=1, then c^(r)=ηc*^(r), where ηis a complex scalar with unit length. It is clear that0≦ρ(c^(r),c*^(r))≦1, and the larger the absolute correlation is, thecloser is the estimated c^(r) to the ground truth. With the absolutecorrelation, the goodness of c′ r=1, 2, . . . , n_(x)n_(y) can bemeasured using the mean absolute correlation (mac) defined as

$\begin{matrix}{{m\; a\; c} = {\sum\limits_{r = 1}^{n_{x}n_{y}}{\frac{\rho \left( {c^{r},c_{*}^{r}} \right)}{n_{x}n_{y}}.}}} & (35)\end{matrix}$

The kernel size is set to k_(x)=k_(y)=6 and the values of the meanabsolute correlation (mac) were explored using different numbers of keptsingular vectors, with the results shown in FIGS. 8( a)-(b). In each ofFIGS. 8( a), (b), and (c), plot 81 represents mac values as calculatedusing a MACO approach according to an embodiment of the invention, plot82 represents mac values as calculated using a B1-Map, and plot 83represents mac values as calculated by dividing the low resolution coilimage by the sum-of-square image. Since A in EQ. (7) is a 288×729matrix, τ can be selected in the interval [1, 288]. As can be observedfrom FIG. 8( a), the value of mac is low, when τ is either very small orvery large. This is reasonable, because on one hand, when τ is verysmall, V_(∥) cannot capture the full range space of the sliding blocks,and thus the left and right hand sides of EQ. (9) do not match well. Onthe other hand, when τ is very large, P=V_(∥)V_(∥) ^(H) defined in EQ.(10) tends to an identity matrix, and in this case, the left and righthand sides of EQ. (9) match very well. However, a near identity matrixfor P gives too much freedom to the sliding blocks in EQ. (12), and thusP contains less information about the calibration data. In the extremecase, when P is an identity matrix, EQ. (12) always holds and it doesnot have any information about the calibration data. Nevertheless, FIGS.8( a)-(b) indicate that one can find a large range of is τ's with whicha MACO approach according to an embodiment of the invention yieldslarger values of mac than both the B1-Map approach and the approachwhich divides the low resolution coil image by the sum-of-square image.The optimal value of τ may be related to the noise in the calibrationdata (note that, for the simulated phantom, no noise is added).Experiments on dozens of data sets show that setting τ as the smallestvalue such that Σ_(k=1) ^(τ)σ_(k)/Σ_(k=1) ²⁸⁸σ_(k)≧0.95 works well,where σ_(k) is the k-th largest singular vector of A in EQ. (7).

In addition, the performance of a MACO approach according to anembodiment of the invention was tested under varying kernel sizes, andthe results are depicted in FIG. 8( c). As can be observed from thefigure, the kernel size can change from 2×2 to 10×10 with k_(x)=k_(y).The value of τ was set to the smallest value for which Σ_(k=1)^(τ)σ_(k)/Σ_(k=1) ²⁸⁸σ_(k)≧0.95. The results again show that a MACOapproach according to an embodiment of the invention can achieve ahigher mac than the B1-Map. The optimal kernel size may depend on thesize of the calibration data, the size of the image, and the noisecontained in the calibration data.

Reconstruction For MR reconstruction, the following formulation wasused:

$\begin{matrix}{{{\min\limits_{m}{\frac{1}{2}{\sum\limits_{i = 1}^{n_{c}}{{{F_{u}\left( {c_{i} \otimes m} \right)} - {yi}}}_{2}^{2}}}} + {\lambda {{Wm}}_{1}}},} & (36)\end{matrix}$

where F_(u) is an under-sampling Fourier transformation operator, c_(i)denotes the coil sensitivity map for the i-th coil, m denotes the imageto be reconstructed, y, is the observed undersampled k-space data forthe i-th coil, and W is a redundant Haar wavelet transformation.

Reconstruction of Two-dimensional Cardiac Image:

Data was acquired from a healthy volunteer on a 1.5T clinical MRscanner. The imaging parameters include: field of view=313×370 mm²,matrix size=176×208, and 20 coils.

FIG. 9 illustrates the reconstruction of a 2D Cardiac Image. Images inthe first row are the coil sensitivity maps estimated by a MACO approachaccording to an embodiment of the invention for coils 1, 10, and 20,respectively; images in the second row are the coil sensitivity mapsestimated by the B1-Map for coils 1, 10, and 20, respectively; the firstimage in the last row illustrates incoherent Cartesian sampling of thek-space; the second image and the third image in the last row correspondto the reconstructions with the CSM's estimated a MACO approachaccording to an embodiment of the invention and the B1-Map,respectively.

The first two rows of FIG. 9 depicts the estimated coil sensitivity mapsby a MACO approach according to an embodiment of the invention andB1-Map with the center 24 lines as calibration data. It can be seen thatthe CSMs estimated by a MACO approach according to an embodiment of theinvention, shown in the first row, are smoother than the B1-Map, shownin the second row. The first image in the last row illustrates theincoherent Cartesian sampling pattern for generating the under-sampledk-space data. Such incoherent sampling leads to an acceleration factorof 6.3. Images obtained by solving EQ. (36) with CSM's estimated by aMACO approach according to an embodiment of the invention, the secondimage in the last row, are smoother than those estimated by the B1-Map,shown in the third image in the last row. The CSM estimated by theB1-Map leads to several artifacts, while the image reconstructed withthe CSM's estimated by a MACO approach according to an embodiment of theinvention contains fewer artifacts due to a more accurate estimation ofthe sensitivity maps.

Reconstruction of Three-Dimensional Abdomen Image:

Data was acquired from a healthy volunteer on a 1.5T clinical MRscanner. Imaging parameters include: field of view=380×320 mm², matrixsize=320×270, 72 partitions, and 24 coils. A Partial Fourier transformis used in the acquisition of the k-space data, and the sampling patternin the phase encoding partition plane is shown in FIG. 10( a). Suchsampling pattern leads to an acceleration factor of 9. An inverseFourier transformation is applied to decouple the samples in the readoutdirection, and the CSM's are estimated and reconstructed for eachdecoupled slice along this direction. Specifically, the CSM's for eachslice are estimated from the center calibration data of size 24×24.Images in the first and second rows of FIG. 10( b) show the coilsensitivity maps estimated by a MACO approach according to an embodimentof the invention and the B1-Map for one partition, and the reconstructedimages for the same partition are shown in the last row of FIG. 10( b).Note that: (1) the CSM's obtained by a MACO approach according to anembodiment of the invention are smoother than those by obtained from aB1-Map; (2) some infolding artifacts are observed in the imagereconstructed with the CSM's estimated by the B1-Map; and (3) the CSM'sestimated by a MACO approach according to an embodiment of the inventionhave fewer artifacts.

Like the B1-Map approach, if c^(r) is the obtained eigenvector by a MACOapproach according to an embodiment of the invention, then ηc^(r) isalso a valid eigenvector for the corresponding eigenvalue if η is acomplex scalar with unit length, i.e., |η|=1. This indicates that a MACOapproach according to an embodiment of the invention may not be able torecover the phase of the reconstructed image. According to embodimentsof the invention, coil sensitivity maps were normalized to the firstcoil, and thus the phase of the reconstructed image is relative to thefirst coil. This may not be an issue when the magnitude of thereconstructed image is the focus.

Influence of the Amount of Calibration Data on the Estimated CSM Quality

In dynamic cardiac MR imaging, the k-space data can be acquired in aninterleaving way, as illustrated in the left plot of FIG. 11, whichillustrates a typical Cartesian sampling pattern used in dynamic cardiacMR imaging. In the figure, solid lines represent acquired data, anddashed lines represent missing data. The average k-space data over time,shown in the right plot of FIG. 11, can be used for calibration purpose.One might think that using more the calibration data will improve thequality of the estimated CSM. However, in real applications, the dataare usually accompanied by noise. More importantly, the boundary highfrequency data might be more sensitive to noise than the center lowfrequency data, as the high frequencies are typically much weaker instrength than the low frequencies, leading to less robustness when thesame noise is added. Therefore, one could conjecture that “the centersquare calibration data could perform better than the center rectangularcalibration data for coil sensitivity maps estimation”.

Experiments were conducted using k-space data acquired on a 1.5Tclinical MR scanner. The acquired k-space data was subsampled togenerate the following under-sampled k-space data: 24 temporal phaseseach containing 18 frequency encodings, 30 coils, and imagesize=192×164. FIGS. 12( a)-(b) illustrate the influence of the amount ofused calibration data on the reconstruction. FIG. 12( a) depicts themask for generating the calibration k-space data and the estimated coilsensitivity maps, where white denotes that the corresponding averagek-space data are used for calibration, and FIG. 12( b) depicts thereconstruction images. FIG. 12( a) illustrates four coil sensitivitymaps estimated by an eigenvector approach according to an embodiment ofthe invention, using (1) the center rectangular 44 lines (top row), (2)the center rectangular 32 lines (2^(nd) row), (3) the center 44×44square block (3^(rd) row), and (4) the center 32×32 square block (bottomrow). It can be seen that, compared with the more rectangularcalibration data, the square center calibration data can achieve asmoother CSM. For example, the black holes 121, 122 in coil 3 and coil 4are observed when using center rectangular calibration data, but notwhen using center square calibration data. FIG. 12( b) compares thereconstruction results achieved by a reconstruction approach accordingto an embodiment of the invention. In FIG. 12( b), the boxed regions 124at the upper left are zoomed out and displayed in the rest of eachimage. The noise observed using a CSM with rectangular calibration datadisappears when using the square calibration data.

In addition, experiments were conducted on synthetic static 2D data asfollows. A phantom of size 256×256 was generated, and then 8 coilsensitivity maps were synthesized from which the k-space data wasgenerated. The CSM was estimated by using (1) the center k_(x)k_(y)square data, and (2) the center rectangular k lines, where k ranges from32 to 64. The mac (maximum correlation) criterion defined above wascomputed. The value of mac is between 0 and 1, and a larger mac valuemeans that the estimated coil profile is closer to the simulated groundtruth. FIG. 13 compares the square calibration data 131 and therectangular calibration data 132 in terms of the mac criterion, fromwhich can be seen a clear advantage of the square calibration data. Theunderlying reason might be that the center low frequency data are lesssensitive to the noise than the boundary high frequency data due to thelarger intensities.

Extension to 3-Dimensions

When the data is acquired using 3D Cartesian sampling, one way forestimating the coil profiles is to decouple the data along thefrequency-encoding direction and then perform the estimation in a 2Dmanner. Then, the 3D coil profiles can be obtained by stacking theestimated 2D coil profiles. The independent estimation of the coilprofiles for each slice causes discontinuity in the profile along thefrequency encoding direction. This discontinuity can be propagated intothe reconstruction, causing artifacts. In addition, applying a 2Deigen-vector approach ignores the correlation between the calibrationdata in the frequency encoding direction and downgrades the effect ofrank reduction in the eigen-vector approach achieved by singular valuedecomposition.

A 3D eigen-vector approach according to an embodiment of the inventionutilizes all the calibration data together and estimates the coilprofiles from different slices simultaneously. FIG. 14 illustrates theconversion of the calibration data into a calibration matrix for 3D CSMestimation, according to an embodiment of the invention. As shown inFIG. 14, the calibration matrix A is constructed from calibration dataof size c_(x)×c_(y)×c_(z) by gathering sliding blocks of kernel sizek_(x)×k_(y)×k_(z). The number of columns in A is[(c_(x)−k_(x)+1)×(c_(y)−k_(y)+1)×(c_(z)−k_(z)+1)], and the number ofrows of A is k_(x)k_(y)k_(z)n_(c), where n_(c) is the number of coils. Alow rank subspace of A can be computer by singular value decompositionA=VΣU^(H), and setting V_(∥) to the left singular vectors of Acorresponding to the leading τ singular values. Let S^(c) be the coilsensitivity at coil cε[1, 2, . . . , c], which are computed following aneigenvector approach according to an embodiment of the invention asdiscussed above, i.e., set to the eigenvectors corresponding to thelargest eigenvalues.

A 3D CSM estimation method according to an embodiment of the inventionwas compared with a 2D method using 3D data. The data were acquired froma healthy volunteer on a 1.5T clinical MR scanner. The imagingparameters include a field of view of 320×260×96 mm², and 26 coils. Inthe phase encoding-partition plane, a total of 2370 (out of 24960)points were acquired, and the central 24×24 block was treated ascalibration data. FIGS. 15( a)-(h) illustrate a comparison of a 2Deigenvector approach and a 3D eigenvector approach according to anembodiment of the invention. FIGS. 15( a), (c), and (e) depict coilprofiles of coils nos. 1, 5 and 10 using a 2D eigenvector approach, andFIGS. 15( b), (d), and (f) depict coil profiles of coils nos. 1, 5 and10 using a 3D eigenvector approach according to an embodiment of theinvention. As shown in FIGS. 15( a)-(f), the vertical direction is thefrequency encoding direction. The coil profiles from the 3D method aresmoother than the ones from the 2D method, especially along thefrequency encoding direction. FIGS. 15( g)-(h) compare the centralregion of the reconstructed image at slice 72 along the transversalaxis, using a 2D eigenvector approach for FIG. 15( g) and a 3Deigenvector approach for FIG. 15( h). The overall image using coilprofiles from a 3D method according to an embodiment of the invention issmoother than the image from the 2D method. The artifact marked by thecircle 152 in FIG. 15( h) is less obvious than the circle 151 in theimage of FIG. 15( g).

Reconstruction with the coil profile from a 3D method according to anembodiment of the invention substantially eliminates artifacts caused bydiscontinuities along the frequency encoding direction from the 2Dmethod. The overall image noise is lower in a 3D method according to anembodiment of the invention due to a smoother coil profile. Efficientcomputation of a 3D eigenvector method according to an embodiment of theinvention can be achieved by utilizing methods according to embodimentsof the invention disclosed above, which has a storage cost ofO(k_(x)n_(y)n_(c) ²).

System Implementations

It is to be understood that embodiments of the present invention can beimplemented in various forms of hardware, software, firmware, specialpurpose processes, or a combination thereof. In one embodiment, thepresent invention can be implemented in software as an applicationprogram tangible embodied on a computer readable program storage device.The application program can be uploaded to, and executed by, a machinecomprising any suitable architecture.

FIG. 16 is a block diagram of an exemplary computer system forimplementing an eigenvector approach to coil sensitivity maps (CSM)estimation for 2-D MR images according to an embodiment of theinvention. Referring now to FIG. 16, a computer system 161 forimplementing the present invention can comprise, inter alia, a centralprocessing unit (CPU) 162, a memory 163 and an input/output (I/O)interface 164. The computer system 161 is generally coupled through theI/O interface 164 to a display 165 and various input devices 166 such asa mouse and a keyboard. The support circuits can include circuits suchas cache, power supplies, clock circuits, and a communication bus. Thememory 163 can include random access memory (RAM), read only memory(ROM), disk drive, tape drive, etc., or a combinations thereof. Thepresent invention can be implemented as a routine 167 that is stored inmemory 163 and executed by the CPU 162 to process the signal from thesignal source 168. As such, the computer system 161 is a general purposecomputer system that becomes a specific purpose computer system whenexecuting the routine 167 of the present invention.

The computer system 161 also includes an operating system and microinstruction code. The various processes and functions described hereincan either be part of the micro instruction code or part of theapplication program (or combination thereof) which is executed via theoperating system. In addition, various other peripheral devices can beconnected to the computer platform such as an additional data storagedevice and a printing device.

It is to be further understood that, because some of the constituentsystem components and method steps depicted in the accompanying figurescan be implemented in software, the actual connections between thesystems components (or the process steps) may differ depending upon themanner in which the present invention is programmed. Given the teachingsof the present invention provided herein, one of ordinary skill in therelated art will be able to contemplate these and similarimplementations or configurations of the present invention.

While the present invention has been described in detail with referenceto exemplary embodiments, those skilled in the art will appreciate thatvarious modifications and substitutions can be made thereto withoutdeparting from the spirit and scope of the invention as set forth in theappended claims.

APPENDIX Proof of Theorem 1:

The proof begins with two technical lemmas.

Lemma 1:

Assume n_(x) and n_(y) are even numbers. Let A, B, P and Q denoten_(x)×n_(y) complex matrices, with the entry at the i-th-row and thej-th column being a_(i,j), b_(i,j), p_(i,j) and q_(i,j), respectively.Let P and Q be the 2-D inverse Fourier transformations of A and B:

$\begin{matrix}{{a_{m,n} = {\sum\limits_{u = 1}^{n_{x}}{\sum\limits_{\upsilon = 1}^{n_{y}}{a_{u,\upsilon}{\exp^{{- 2}\pi \; {j{({{{({u - 1})}{m/n_{x}}} + {{({\upsilon - 1})}{n/n_{y}}}})}}}/\sqrt{n_{x}n_{y}}}}}}},} & (37) \\{{q_{m,n} = {\sum\limits_{u = 1}^{n_{x}}{\sum\limits_{\upsilon = 1}^{n_{y}}{b_{u,\upsilon}{\exp^{{- 2}\pi \; {j{({{{({u - 1})}{m/n_{x}}} + {{({\upsilon - 1})}{n/n_{y}}}})}}}/\sqrt{n_{x}n_{y}}}}}}},{{{where}\mspace{14mu} j} = {\sqrt{- 1}.\mspace{14mu} {If}}}} & (38) \\{{{a_{u,\upsilon} = \overset{\_}{b_{{n_{x} + 2 - u},{n_{y} + 2 - \upsilon}}}},{u = 2},\ldots \mspace{14mu},n_{x},{\upsilon = 2},\ldots \mspace{14mu},n_{y}}{and}} & (39) \\{{q_{1,\upsilon} = {a_{u,1} = {b_{1,\upsilon} = {b_{u,1} = 0}}}},{u = 1},\ldots \mspace{14mu},n_{x},{\upsilon = 1},\ldots \mspace{14mu},n_{y},{{{one}\mspace{14mu} {has}\mspace{14mu} p_{m,n}} = \overset{\_}{q_{m,n}}},{m = 1},\ldots \mspace{14mu},n_{x},{n = 1},\ldots \mspace{14mu},{n_{y}.}} & (40)\end{matrix}$

Proof:

Inserting EQS. (39) and (40) into EQ. (37), one obtains EQ. (41):

$\begin{matrix}\begin{matrix}{p_{m,n} = {\sum\limits_{u = 2}^{n_{x}}{\sum\limits_{\upsilon = 2}^{n_{y}}{\overset{\_}{b_{{n_{x} + 2 - u},{n_{y} + 2 - \upsilon}}}{\exp^{{- 2}\pi \; {j{({{{({u - 1})}{m/n_{x}}} + {{({\upsilon - 1})}{n/n_{y}}}})}}}/\sqrt{n_{x}n_{y}}}}}}} \\{= {\sum\limits_{u^{\prime} = 2}^{n_{x}}{\sum\limits_{\upsilon^{\prime} = 2}^{n_{y}}{\overset{\_}{b_{u^{\prime},\upsilon^{\prime}}}{\exp^{{- 2}\pi \; {j{({{{({n_{x} + 1 - u^{\prime}})}{m/n_{x}}} + {{({n_{y} + 1 - \upsilon^{\prime}})}{n/n_{y}}}})}}}/\sqrt{n_{x}n_{y}}}}}}} \\{= {\sum\limits_{u^{\prime} = 2}^{n_{x}}{\sum\limits_{\upsilon^{\prime} = 2}^{n_{y}}{\overset{\_}{b_{u^{\prime},\upsilon^{\prime}}}{\exp^{{- 2}\pi \; {j{({{{({1 - u^{\prime}})}{m/n_{x}}} + {{({1 - \upsilon^{\prime}})}{n/n_{y}}}})}}}/\sqrt{n_{x}n_{y}}}}}}} \\{= {\sum\limits_{u^{\prime} = 2}^{n_{x}}{\sum\limits_{\upsilon^{\prime} = 2}^{n_{y}}{\overset{\_}{b_{u^{\prime},\upsilon^{\prime}}\exp^{{- 2}\pi \; {j{({{{({u^{\prime} - 1})}{m/n_{x}}} + {{({\upsilon^{\prime} - 1})}{n/n_{y}}}})}}}}/\sqrt{n_{x}n_{y}}}}}} \\{= {\sum\limits_{u^{\prime} = 1}^{n_{x}}{\sum\limits_{\upsilon^{\prime} = 1}^{n_{y}}{\overset{\_}{b_{u^{\prime},\upsilon^{\prime}}\exp^{{- 2}\pi \; {j{({{{({u^{\prime} - 1})}{m/n_{x}}} + {{({\upsilon^{\prime} - 1})}{n/n_{y}}}})}}}}/\sqrt{n_{x}n_{y}}}}}} \\{= {\overset{\_}{q_{m,n}}.}}\end{matrix} & (41)\end{matrix}$

Lemma 2:

Let u and v denote k_(x)k_(y)×1 column vectors with complex entries. Letu_(t) and v_(t) denote the t-th elements of u and v, respectively.Define

$\begin{matrix}{{p = {\sum\limits_{t = 1}^{k_{x}k_{y}}{\mathcal{F}^{H}\left( {_{t}\left( {\overset{\_}{u_{t}}v} \right)} \right)}}},} & (42) \\{{q = {\sum\limits_{t = 1}^{k_{x}k_{y}}{\mathcal{F}^{H}\left( {_{t}\left( {\overset{\_}{\upsilon_{t}}u} \right)} \right)}}},} & (43)\end{matrix}$

which are n_(x)n_(y)×1 column vectors. Then,

p= q.  (44)

Proof:

Let

$\begin{matrix}{{a = {\sum\limits_{k = 1}^{k_{x}k_{y}}{_{t}\left( {\overset{\_}{u_{t}}v} \right)}}},} & (45) \\{b = {\sum\limits_{t = 1}^{k_{x}k_{y}}{{_{t}\left( {\overset{\_}{\upsilon_{t}}u} \right)}.}}} & (46)\end{matrix}$

Let A and B be n_(x)×n_(y) complex matrices that are the matrix versionof a and b respectively. Making use of the definition of the zeropadding operator, and

a _(u,v) =b _(n) _(x) _(+2−u,n) _(y) _(+2−v) ,u=2, . . . , n _(x) ,v=2,. . . , n _(y),  (47)

and

a _(1,v) =a _(u,1) =b _(1,v) =b _(u,1)=0,u=0,u=1, . . . , n _(x) ,v=1, .. . , n _(y).  (48)

Applying Lemma 2, one can verify EQ. (44).

Theorem 1 can now be proved. Denote q_(j,i) ^(r) as the j-th row and thei-th column of S^(r)M^(H). It follows from EQS. (17) and (18) that

$\begin{matrix}{{q_{j,i}^{r} = {\sum\limits_{t = 1}^{k_{x}k_{y}}s_{i,j,t}^{r}}},{\forall{r.}}} & (49)\end{matrix}$

Denote q_(j,i) as an n_(x)×n_(y) column vector with the r-th entry beingq_(j,i) ^(r). Then

$\begin{matrix}{q_{j,i} = {\sum\limits_{t = 1}^{k_{x}k_{y}}{s_{i,j,t}.}}} & (50)\end{matrix}$

To prove that S^(r)M^(H) is Hermitian, it suffices to verify

q _(j,i)= q _(i,j)   (51)

It follows from EQS. (9) and (10) that

$\begin{matrix}{{p_{i,j,t} = {\sum\limits_{k = 1}^{\tau}{\overset{\_}{\upsilon_{i,k}^{t}}v_{j,k}}}},} & (52)\end{matrix}$

where v_(i,k) ^(t) denotes the t-th element in v_(i,k). Therefore,

$\begin{matrix}\begin{matrix}{q_{j,i} = {\sum\limits_{t = 1}^{k_{x}k_{y}}{\mathcal{F}^{H}\left( {_{t}\left( p_{i,j,t} \right)} \right)}}} \\{= {\sum\limits_{t = 1}^{k_{x}k_{y}}{\mathcal{F}^{H}\left( {_{t}\left( {\sum\limits_{k = 1}^{\tau}{\overset{\_}{\upsilon_{i,k}^{t}}v_{j,k}}} \right)} \right.}}} \\{= {\sum\limits_{k = 1}^{\tau}\left( {\sum\limits_{t = 1}^{k_{x}k_{y}}{{\mathcal{F}^{H}\left( {_{t}\left( {\overset{\_}{\upsilon_{i,k}^{t}}v_{j,k}} \right)} \right)}.}} \right.}}\end{matrix} & (53)\end{matrix}$

By utilizing Lemma 2, one can obtain EQ (51).

Since

$\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}$

is Hermitian, its eigenvalues are real. This ends the proof of thistheorem.

Proof of Theorem 2:

For any pair of λ and α that satisfy EQ. (20), left multiplying bothsides of EQ. (20) by M yields

M(S ^(r))^(H) α=λMM ^(H) α=λk _(x) k _(y)α  (54)

where MM^(H)=k_(x)k_(y)I was used for establishing the last equality.Therefore, λ and α are also the eigenvalue and eigenvector of EQ. (21),respectively.

Proof of Theorem 3:

For any pair of λ and α that satisfy EQ. (20), left multiplying bothsides of EQ. (20) by S^(r) yields

S ^(r)(S ^(r))^(H) α=λS ^(r) M ^(H)α=λ¹ k _(x) k _(y)α  (55)

where Theorem 2 establishes the second eqaulity. Therefore, λ and α arealso the eigenvalue and eigenvector of EQ. (22), respectively.

For any pair of λ and α that satisfy EQ. (25),

S ^(r)(S ^(r))^(H) MM ^(H)α=(λk _(x) k _(y))² α=S ^(r) M ^(H) S ^(r) M^(H) α  (56)

where MM^(H)=k_(x)k_(y)I and (S^(r))^(H)M=M^(H)S^(r) have been used (seeTheorem 1). The last equality in EQ. (56) indicates that α is theeigenvector of S^(r)M^(H)S^(r)M^(H) corresponding to eigenvalue(λk_(x)k_(y))².

Proof of Theorem 4:

The proof begins with a technical lemma, whose proof is omitted here.

Lemma 3: F^(H)(P_(t)(v_(j,k)))

F^(H)(P_(t)(v_(j′,k′))) is invariant to the spatial location t, i.e.,F^(H)(P_(t)(v_(j,k)))

F^(H)(P_(t)(v_(j′,k′)))=F^(H)(P_(t′)(v_(j,k)))

F^(H)(P_(t′)(v_(j′,k′))).

Now, the theorem may be proved. The j-th row and j′-th column ofS^(r)(S^(r))^(H) can be computed as

$\begin{matrix}{{\sum\limits_{i = 1}^{n_{c}}m_{i,j,j^{\prime}}^{r}},{where}} & (57) \\{{m_{i,j,j^{\prime}}^{r} = {\sum\limits_{t = 1}^{k_{x}k_{y}}{s_{i,j,t}^{r}\overset{\rightharpoonup}{s_{i,j^{\prime},t}^{r}}}}},{\forall{r.}}} & (58)\end{matrix}$

EQ. (58) can be rewritten as the following matrix format:

$\begin{matrix}{m_{i,j,j^{\prime}} = {\sum\limits_{t = 1}^{k_{x}k_{y}}{s_{i,j,t} \otimes {\overset{\_}{s_{i,j^{\prime},t}}.}}}} & (59)\end{matrix}$

Inserting EQ. (52) into EQ. (15) yields

$\begin{matrix}{s_{i,j,t} = {{F^{H}\left( {P_{t}\left( p_{i,j,t} \right)} \right)} = {\sum\limits_{k = 1}^{\tau}{\overset{\_}{v_{i,k}^{t}}{F^{H}\left( {P_{t}\left( v_{j,k} \right)} \right)}}}}} & (60)\end{matrix}$

In computing s_(i,j,t)

s_(i,j′,t) , a useful operation is

F ^(H)(P _(t)(v _(j,k)))

F ^(H)(P _(t)(v _(j′,k′)))  (61)

which is identical to

F ^(H)(P _(c)(v _(j,k)))

F^(H)(P _(c)(v _(j′,k′)))  (62)

using Lemma 3.Let Γ_(j,j′) ^(r) be a τ×τ matrix whose k-th row and k′-th column arethe value of z_(j,k) ^(r) z_(j′,k′) ^(r) . Denote v_(i,:) as the i-throw of V_(∥). Then

m _(i,j,j′) ^(r) =k _(x) k _(y) v _(i,:) Γ_(j,j′) ^(r) v _(i,:) ^(H).  (63)

It is interesting to note that

$\begin{matrix}\begin{matrix}{{\sum\limits_{i = 1}^{n_{c}}m_{i,j,j^{\prime}}^{r}} = {k_{x}k_{y}{\sum\limits_{i = 1}^{n_{c}}{\overset{\_}{v_{i,:}}\Gamma_{j,j^{\prime}}^{r}\overset{\_}{v_{i,:}^{H}}}}}} \\{= {k_{x}k_{y}{\sum\limits_{i = 1}^{n_{c}}{{trace}\left( {\Gamma_{j,j^{\prime}}^{r}\overset{\_}{v_{i,:}^{H}v_{i,:}}} \right)}}}} \\{= {k_{x}k_{y}{{trace}\left( {\Gamma_{j,j^{\prime}}^{r}{\sum\limits_{i = 1}^{n_{c}}\overset{\_}{v_{i,:}^{H}v_{i,:}}}} \right)}}} \\{= {k_{x}k_{y}{{trace}\left( \Gamma_{j,j^{\prime}}^{r} \right)}}} \\{{= {k_{x}k_{y}{\sum\limits_{k = 1}^{\tau}{z_{j,k}^{r}\overset{\_}{z_{j,k}^{r}}}}}},}\end{matrix} & (64)\end{matrix}$

where trace(A) denotes the summation of the diagonal entries of a matrixA. Therefore, EQ. (32) holds.

What is claimed is:
 1. A method for estimating a coil sensitivity mapfor a magnetic resonance (MR) image, comprising the steps of: providinga matrix A of sliding blocks of a 2D n_(x)×n_(y) image of coilcalibration data, wherein ${A = \begin{pmatrix}a_{1,1} & a_{1,2} & \ldots & a_{1,n_{b}} \\a_{2,1} & a_{2,2} & \ldots & a_{2,n_{b}} \\\ldots & \ldots & \; & \ldots \\a_{n_{c},1} & a_{n_{c},2} & \ldots & a_{n_{c},n_{b}}\end{pmatrix}},$ n_(c) is a number of coils, n_(b) is a number ofsliding blocks extracted from the coil calibration data, and a_(i,j) isa k_(x)k_(y)×1 column vector that represents a jth sliding block of anith coil; determining a unit eigenvector α that maximizes

(S^(r))^(H)α,M^(H)α

, wherein S^(r) is an inverse Fourier transform of a zero-padded matrixP=V_(∥)V_(∥) ^(H) at spatial location r in the 2D n_(x)×n_(y) image,wherein$V_{} = {\left\lbrack {v_{1\;},v_{2},\ldots \mspace{14mu},v_{\tau}} \right\rbrack = \begin{pmatrix}v_{1,1} & v_{1,2} & \ldots & v_{1,\tau} \\v_{2,1} & v_{2,2} & \ldots & v_{2,\tau} \\\ldots & \ldots & \; & \ldots \\v_{n_{c},1} & v_{n_{c},2} & \ldots & v_{n_{c},\tau}\end{pmatrix}}$ is a matrix composed of left singular vectors of Acorresponding to τ leading singular values wherein v_(i,k) is ak_(x)k_(y)×1 column vector that is an i-th block in vector v_(k),$M = \left( {\begin{pmatrix}1 & 1 & \ldots & 1 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\begin{pmatrix}0 & 0 & \ldots & 0 \\1 & 1 & \ldots & 1 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\1 & 1 & \ldots & 1\end{pmatrix}} \right)$ is a n_(c)×k_(x)k_(y)n_(c) sparse matrix of thesame size as S^(r); and finding an optimizer c^(r) that maximizes acorrelation between (S^(r))^(H)α and M^(H)α wherein said optimizer c^(r)is a column vector of coil sensitivities of a coil at spatial position rin the 2D n_(x)×n_(y) image.
 2. A method for estimating a coilsensitivity map for a magnetic resonance (MR) image, comprising thesteps of: providing a matrix A of sliding blocks of a 2D n_(x)×n_(y)image of coil calibration data, wherein ${A = \begin{pmatrix}a_{1,1} & a_{1,2} & \ldots & a_{1,n_{b}} \\a_{2,1} & a_{2,2} & \ldots & a_{2,n_{b}} \\\ldots & \ldots & \; & \ldots \\a_{n_{c},1} & a_{n_{c},2} & \ldots & a_{n_{c},n_{b}}\end{pmatrix}},$ n_(c) is a number of coils, n_(b) is a number ofsliding blocks extracted from the coil calibration data, and a_(i,j) isa k_(x)k_(y)×1 column vector that represents a jth sliding block of anith coil; calculating a left singular matrix V_(∥) from a singular valuedecomposition of A, wherein A=VΣU^(H) and$V_{} = {\left\lbrack {v_{1},v_{2},\ldots \mspace{14mu},v_{\tau}} \right\rbrack = \begin{pmatrix}v_{1,1} & v_{1,2} & \ldots & v_{1,\tau} \\v_{2,1} & v_{2,2} & \ldots & v_{2,\tau} \\\ldots & \ldots & \; & \ldots \\v_{n_{c},1} & v_{n_{c},2} & \ldots & v_{n_{c},\tau}\end{pmatrix}}$ is a matrix of left singular vectors of A correspondingto τ leading singular values wherein v_(i,k) is a k_(x)k_(y)×1 columnvector that is an i-th block in vector v_(k); calculating$\begin{matrix}{P = {V_{}V_{}^{H}}} \\{= \left( {\begin{pmatrix}p_{1,1,1} & p_{1,1,2} & \ldots & p_{1,1,{k_{x}k_{y}}} \\p_{1,1,1} & p_{1,2,2} & \ldots & p_{1,2,{k_{x}k_{y}}} \\\ldots & \ldots & \ldots & \ldots \\p_{1,n_{c},1} & p_{1,n_{c},2} & \ldots & p_{1,n_{c},{k_{x}k_{y}}}\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}p_{n_{c},1,1} & \ldots & p_{n_{c},1,{k_{x}k_{y}}} \\p_{n_{c},2,1} & \ldots & p_{n_{c},2,k_{x},k_{y}} \\\ldots & \ldots & \ldots \\p_{n_{c},n_{c},1} & \ldots & p_{n_{c},n_{c},k_{x},k_{y}}\end{pmatrix}} \right)}\end{matrix}$ wherein p_(i,j,t) is a k_(x)k_(y)×1 column vector;calculating S_(i,j,t)=F^(H)(P_(t)(p_(i,j,t))) wherein F^(H) representsan inverse Fourier transform and P_(t) represents a zero-paddingoperator; and solving M^(H)c^(r) (S^(r))^(H)c^(r) for c^(r), where c^(r)is a vector of coil sensitivity maps for all coils at spatial locationr, $S^{r} = \left( {\begin{pmatrix}s_{1,1,1}^{r} & s_{1,1,2}^{r} & \ldots & s_{1,1,{k_{x}k_{y}}}^{r} \\s_{1,2,1}^{r} & s_{1,2,2}^{r} & \; & s_{1,2,{k_{x}k_{y}}}^{r} \\\ldots & \ldots & \ldots & \ldots \\s_{1,n_{c},1}^{r} & s_{1,n_{c},2}^{r} & \ldots & s_{1,n_{c},{k_{x}k_{y}}}^{r}\end{pmatrix}\mspace{20mu} \ldots \mspace{20mu} \begin{pmatrix}s_{n_{c},1,1}^{r} & \ldots & s_{n_{c},1,{k_{x}k_{y}}}^{r} \\s_{n_{c},2,1}^{r} & \ldots & s_{n_{c},2,k_{x},k_{y}}^{r} \\\ldots & \ldots & \ldots \\s_{n_{c},n_{c},1}^{r} & \ldots & s_{n_{c},n_{c},k_{x},k_{y}}^{r}\end{pmatrix}} \right)$   and $M = {\left( {\begin{pmatrix}1 & 1 & \ldots & 1 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\begin{pmatrix}0 & 0 & \ldots & 0 \\1 & 1 & \ldots & 1 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\1 & 1 & \ldots & 1\end{pmatrix}} \right).}$
 3. The method of claim 2, further comprisingfinding an optimizer c_(r) that maximizes a correlation between(S^(r))^(H)α and M^(H)α wherein said optimizer c^(r) is vector of coilsensitivity maps for all coils at spatial location r in the 2Dn_(x)×n_(y) image.
 4. The method of claim 2, wherein solvingM^(H)c^(r)=(S_(r))^(H)c^(r) comprises solving eigenvalue equations${{\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}\beta} = {\overset{\_}{\lambda}\; \beta}},{{{{and}\left( {S^{r}\left( S^{r} \right)}^{H} \right)}\gamma} = {\overset{\sim}{\lambda \;}\gamma}},$wherein λ and {tilde over (λ)} denote eigenvalues, and β and γ denoteeigenvectors.
 5. The method of claim 4, further comprising calculatingS^(r)M^(H) using${m_{i,j} = {{\sum\limits_{t = 1}^{k_{x}k_{y}}s_{i,j,t}} = {{\sum\limits_{t = 1}^{k_{x}k_{y}}{F^{H}\left( {P_{t}\left( p_{i,j,t} \right)} \right)}} = {F^{H}\left( {\sum\limits_{t = 1}^{k_{x}k_{y}}{P_{t}\left( p_{i,j,t} \right)}} \right)}}}},$wherein m_(i,j) is an n_(x)n_(y)×1 column vector containing the entriesof m_(i,j) ^(r) at all spatial locations, and m_(i,j) ^(r) is the i,jelements of the 2D MR image at spatial location r.
 6. The method ofclaim 4, further comprising calculating S^(r)(S^(r))^(H) asS^(r)(S_(r))^(H)=k_(x)k_(y)Z^(r)(Z^(r))^(H), wherein$Z^{r} = \begin{pmatrix}z_{1,1}^{r} & z_{1,2}^{r} & \; & z_{1,\tau}^{r} \\z_{2,1}^{r} & z_{2,2}^{r} & \; & z_{2,\tau}^{r} \\\; & \; & \; & \; \\z_{n_{c},1}^{r} & z_{n_{c},2}^{r} & \; & z_{n_{c},\tau}^{r}\end{pmatrix}$ is an n_(c)×τ matrix with the j-row and k-th column beingz_(j,k) ^(r), z_(j,k) ^(r) denotes the r-th entry ofz_(j,k)=F^(H)(P_(c)(v_(j,k))), and P_(c) denotes the zero-paddingoperator that maps a center of v_(j,k) to a center of z_(j,k).
 7. Themethod of claim 6, further comprising computing the i-th row and j-thcolumn of${Z^{r}\left( Z^{r} \right)}^{H}\mspace{14mu} {as}\mspace{14mu} {\sum\limits_{k = 1}^{\tau}{z_{i,k}^{r}\overset{\_}{z_{i,k}^{r}}}}$wherein z_(i,k) ^(r) denotes the i-th row and k-th column in Z_(r). 8.The method of claim 2, wherein the coil calibration data is acquiredfrom a center square of frequency space data.
 9. The method of claim 2,wherein the matrix A is constructed from calibration data of sizec_(x)×c_(y)×c_(z) by gathering sliding blocks of kernel sizek_(x)×k_(y)×k_(z), wherein A has[(c_(x)−k_(x)+1)×(c_(y)−k_(y)+1)×(c_(z)−k_(z)+1)] columns andk_(x)k_(y)k_(z)n_(c) rows.
 10. A non-transitory program storage devicereadable by a computer, tangibly embodying a program of instructionsexecuted by the computer to perform the method steps for estimating acoil sensitivity map for a magnetic resonance (MR) image, the methodcomprising the steps of: providing a matrix A of sliding blocks of a 2Dn_(x)×n_(y) image of coil calibration data, wherein${A = \begin{pmatrix}a_{1,1} & a_{1,2} & \ldots & a_{1,n_{b}} \\a_{2,1} & a_{2,2} & \ldots & a_{2,n_{b}} \\\ldots & \ldots & \; & \ldots \\a_{n_{c},1} & a_{n_{c},2} & \ldots & a_{n_{c},n_{b}}\end{pmatrix}},$ n_(c) is a number of coils, n_(b) is a number ofsliding blocks extracted from the coil calibration data, and a_(i,j) isa k_(x)k_(y)×1 column vector that represents a jth sliding block of anith coil; calculating a left singular matrix V_(∥) from a singular valuedecomposition of A, wherein A=VΣU^(H) and$V_{} = {\left\lbrack {v_{1},v_{2},\ldots \mspace{14mu},v_{\tau}} \right\rbrack = \begin{pmatrix}v_{1,1} & v_{1,2} & \ldots & v_{1,\tau} \\v_{2,1} & v_{2,2} & \ldots & v_{2,\tau} \\\ldots & \ldots & \; & \ldots \\v_{n_{c},1} & v_{n_{c},2} & \ldots & v_{n_{c},\tau}\end{pmatrix}}$ is a matrix of left singular vectors of A correspondingto τ leading singular values wherein v_(i,k) is a k_(x)k_(y)×1 columnvector that is an i-th block in vector v_(k); calculating$\begin{matrix}{P = {V_{}V_{}^{H}}} \\{= \left( {\begin{pmatrix}p_{1,1,1} & p_{1,1,2} & \ldots & p_{1,1,{k_{x}k_{y}}} \\p_{1,1,1} & p_{1,2,2} & \ldots & p_{1,2,{k_{x}k_{y}}} \\\ldots & \ldots & \ldots & \ldots \\p_{1,n_{c},1} & p_{1,n_{c},2} & \ldots & p_{1,n_{c},{k_{x}k_{y}}}\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}p_{n_{c},1,1} & \ldots & p_{n_{c},1,{k_{x}k_{y}}} \\p_{n_{c},2,1} & \ldots & p_{n_{c},2,k_{x},k_{y}} \\\ldots & \ldots & \ldots \\p_{n_{c},n_{c},1} & \ldots & p_{n_{c},n_{c},k_{x},k_{y}}\end{pmatrix}} \right)}\end{matrix}$ wherein p_(i,j,t) is a k_(x)k_(y)×1 column vector;calculating S_(i,j,t)=F^(H)(P_(t)(p_(i,j,t))) wherein F^(H) representsan inverse Fourier transform and P_(t) represents a zero-paddingoperator; and solving M^(H)c^(r)=(S^(r))^(H)c_(r) for c^(r), where c^(r)is a vector of coil sensitivity maps for all coils at spatial locationr, $S^{r} = \left( {\begin{pmatrix}s_{1,1,1}^{r} & s_{1,1,2}^{r} & \ldots & s_{1,1,{k_{x}k_{y}}}^{r} \\s_{1,2,1}^{r} & s_{1,2,2}^{r} & \; & s_{1,2,{k_{x}k_{y}}}^{r} \\\ldots & \ldots & \ldots & \ldots \\s_{1,n_{c},1}^{r} & s_{1,n_{c},2}^{r} & \ldots & s_{1,n_{c},{k_{x}k_{y}}}^{r}\end{pmatrix}\mspace{20mu} \ldots \mspace{20mu} \begin{pmatrix}s_{n_{c},1,1}^{r} & \ldots & s_{n_{c},1,{k_{x}k_{y}}}^{r} \\s_{n_{c},2,1}^{r} & \ldots & s_{n_{c},2,k_{x},k_{y}}^{r} \\\ldots & \ldots & \ldots \\s_{n_{c},n_{c},1}^{r} & \ldots & s_{n_{c},n_{c},k_{x},k_{y}}^{r}\end{pmatrix}} \right)$   and $M = {\left( {\begin{pmatrix}1 & 1 & \ldots & 1 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\begin{pmatrix}0 & 0 & \ldots & 0 \\1 & 1 & \ldots & 1 \\\ldots & \ldots & \; & \ldots \\0 & 0 & \ldots & 0\end{pmatrix}\mspace{14mu} \ldots \mspace{14mu} \begin{pmatrix}0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 \\\ldots & \ldots & \; & \ldots \\1 & 1 & \ldots & 1\end{pmatrix}} \right).}$
 11. The computer readable program storagedevice of claim 10, the method further comprising finding an optimizerc^(r) that maximizes a correlation between (S^(r))^(H)α and M^(H)αwherein said optimizer c^(r) is vector of coil sensitivity maps for allcoils at spatial location r in the 2D n_(x)×n_(y) image.
 12. Thecomputer readable program storage device of claim 10, wherein solvingM^(H)c^(r)=(S^(r))^(H)c^(r) comprises solving eigenvalue equations${{\frac{{M\left( S^{r} \right)}^{H}}{k_{x}k_{y}}\beta} = {\overset{\_}{\lambda}\; \beta}},$and(S ^(r)(S ^(r))^(H))γ={tilde over (λ)}γ, wheren λ and {tilde over (λ)}denote eigenvalues, and β and γ denote eigenvectors.
 13. The computerreadable program storage device of claim 12, the method furthercomprising calculating S^(r)M^(H) using${m_{i,j} = {{\sum\limits_{t = 1}^{k_{x}k_{y}}s_{i,j,t}} = {{\sum\limits_{t = 1}^{k_{x}k_{y}}{F^{H}\left( {P_{t}\left( p_{i,j,t} \right)} \right)}} = {F^{H}\left( {\sum\limits_{t = 1}^{k_{x}k_{y}}{P_{t}\left( p_{i,j,t} \right)}} \right)}}}},$wherein m_(i,j) is an n_(x)n_(y)×1 column vector containing the entriesof m_(i,j) ^(r) at all spatial locations, and m_(i,j) ^(r) is the i,jelements of the 2D MR image at spatial location r.
 14. The computerreadable program storage device of claim 12, the method furthercomprising calculating S^(r)(S^(r))^(H) asS^(r)(S^(r))^(H)=k_(x)k_(y)Z^(r)(Z^(r))^(H), wherein$Z^{r} = \begin{pmatrix}z_{1,1}^{r} & z_{1,2}^{r} & \ldots & z_{1,\tau}^{r} \\z_{2,1}^{r} & z_{2,2}^{r} & \ldots & z_{2,\tau}^{r} \\\ldots & \ldots & \; & \ldots \\z_{n_{c},1}^{r} & z_{n_{c},2}^{r} & \ldots & z_{n_{c},\tau}^{r}\end{pmatrix}$ is an n_(c)×τ matrix with the j-row and k-th column beingz_(j,k) ^(r), z_(j,k) ^(r) denotes the r-th entry ofz_(j,k)=F^(H)(P^(c)(v_(j,k))), and P_(c) denotes the zero-paddingoperator that maps a center of v_(j,k) to a center of z_(j,k).
 15. Thecomputer readable program storage device of claim 14, the method furthercomprising computing the i-th row and j-th column of Z^(r)(Z^(r))^(H) as$\sum\limits_{k = 1}^{\tau}{z_{i,k}^{r}\overset{\_}{z_{i,k}^{r}}}$wherein z_(i,k) ^(r) denotes the i-th row and k-th column in Z^(r). 16.The computer readable program storage device of claim 10, wherein thecoil calibration data is acquired from a center square of frequencyspace data.
 17. The computer readable program storage device of claim10, wherein the matrix A is constructed from calibration data of sizec_(x)×c_(y)×c_(z) by gathering sliding blocks of kernel sizek_(x)×k_(y)×k_(z), wherein A has[(c_(x)−k_(x)+1)×(c_(y)−k_(y)+1)×(c_(z)−k_(z)+1)] columns andk_(x)k_(y)k_(z)n_(c) rows.